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The Ehipsoidal Statistical model (ES-model) and the Shakhov model (S-model) are con- 
structed for the correction of Prandtl number of the original BGK model through the modi- 



' O , fication of stress and heat flux. Even though in the continuum flow regime, both models can 

^ , give the same Navier-Stokes equations with correct Prandtl number, their modiflcation of 

• ' the collision term may have different dynamic effect in the non-equilibrium transition flow 

. ^z^ ' regimes. With the introduction of one free parameter, a generalized kinetic model with the 

^-^, combination of the ES-model and S-model can be developed, and this new model can get 

Oj' the correct Navier-Stokes equations in the continuum flow regime as well, but with abundant 

dynamic effect through the adjustment of the new degree of freedom. In order to validate 

<^ ' the generalized model, a numerical method based on the unifled gas kinetic scheme (UGKS) 

/^ I has been developed for the new model. The physical performance of the new model with 

00 , the variation of the free parameter has been tested, where the ES-model and S-model be- 

l' ' come the limiting cases. In transition flow regime, many physical problems, i.e., the shock 

^— ^ ' structure and micro-flows, have been studied using the generalized model. With a careful 

choice of the free parameter, good results can be achieved for most test cases. The overall 



conclusion is that the S-model predicts more accurate numerical solutions in most tough 
test cases presented in this paper than the ES-model, while ES-model performs better in the 



X 

.\^.' cases when the flow is mostly driven by heat, such as a channel flow with large boundary 



temperature variations at high Knudsen number. The numerical study demonstrates the 
necessity of developing such a generalized model. With the inclusion of one more freedom, 
in the transition regime the new kinetic model may provide more accurate solution than the 
ES and Shakhov models. 
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I. INTRODUCTION 

The monatomic rarefied gas behavior can be described by the Boltzmann equation. However, 
the colhsion term of the Boltzmann equation is a multiple integral term which is very complicated 
for analysis and numerical computation. The kinetic rnodel is a simplification of the Boltzmann 
equation. The simplest kinetic model is the BGK model [l!] in which the collision term is replaced by 
a relaxation term. This relaxation term in the BGK model mimics the main relaxation process from 
nonequilibrium state towards to a local equilibrium one with a Maxwellian distribution function. 
The local equilibrium state is determined by the local conservative flow variables, namely, the 
density, the momentum and the energy. The BGK model becomes an important kinetic model 
for analysis and numerical simulation of nonequilibrium flows. However, the Chapman-Enskog 
expansion of the BGK model gives the Navier-Stokes equations with a unit Prandtl number, which 
is different from the physical reality in the continuum flow regime. For a monatomic gas, the 
accepted Prandtl number is about 2/3 in a wide range of flow conditions. 

In order to flx the Prandtl number, many kinetic models have been proposed in the past 
decades. The main idea is to modify the relaxation states. For example, the Ellipsoidal Statistical 
BGK model [2] employs a Gaussian distribution as the relaxation equilibrium state instead of the 
Maxwellian. This model is not very popular until the proof of the entropy condition by Andries 
et. al. [3|. In the ES-model, besides the conservative flow variables, the local stress tensor also 
involves in the post-collision state. By changing the free parameter in the ES-model, it can present 
an arbitrary Prandtl number. Moreover, the nonnegative property of the Gaussian distribution 
becomes a favorable physical property. 

Another very popular kinetic model is the Shakhov model [4|. Unlike the ES-model, it adjusts 
the heat flux in the relaxation term, but keeps the stress tensor the same as the original BGK one. 
The Hermit polynomial is adopted to modify the heat flux. So in terms of low order moments, the 
S-model keeps the same as the BGK model. The S-model also presents a correct Prandtl number. 
But, it allows negative value of distribution function, and its H-theorem was only proved in near 



i 



equilibrium condition 

In 1990, Liu proposed a new kinetic model by considering the gain term and lost term of the 
Boltzmann equation separately [5]. He used the Chapman-Enskog distribution directly to evaluate 
the relaxation term. The modification of the collision term involves the space derivatives. Liu 
model changes both the heat flux and stress tensor of the relaxation process, and provides a correct 
Prandtl number in the continuum flow regime. Due to its relatively complicated formulation, this 



model has not been widely used. 

Although all above models provide correct Prandtl number in the continuum flow regime, their 
properties are very different in the transition regime |6Hl0l|. Garzo reported a singular behavior 
of Liu model and attributed it to the negative distribution function [£]. Graur studied the heat 
transfer problem, and found that the ES-model provides better results than the S-model through 
the comparison with the results from the Boltzmann equation. The ES-model keeps the distribution 
function positive, while the S-model and Liu model always allow un-physical negative distribution 
function. It seems that the nonnegative properties of the ES-model are important and promising. 
Moreover, the ES-model satisfies the H theorem, while the H-theorem of the S-model is only proved 
in the near local equilibrium state [7|]. ^However, some other studies did not tell the same story. 



n 

Mieussens [Ul] and Kudryavtsev et al. [6] both reported the early rising of temperature profile in 



the shock structure solution by the ES-model. 

In fact, the physical performance of these models has not yet been evaluated extensively in 
the transition regime. The properties, such as the H-theorem, nonnegative distribution, and con- 
servation etc., cannot cover a complete picture of dynamics of the particle collision term and the 
evolution of the distribution function. Since the original motivation for the development of the 
kinetic models is to fix the Prandtl number which is well defined in the continuum fiow regime, in 
transition regime it is expected that significant differences in their performance would appear in dif- 
ferent physical problems. Furthermore, the practical requirement cares more about the macroscopic 
quantities, such as the moments of a distribution function. The H-theorem and the nonnegative 
distribution function cannot guarantee a correct dynamic evolution of macroscopic quantities. So, 
it is necessary to inspect the practical performance of different kinetic models through the numerical 
simulations in the transition regime. In order to cover a whole spectrum of dynamic performance of 
kinetic models, we are going to introduce a generalized kinetic model which combines the ES-model 
and S-model. With the combination of these two models, besides the correct capturing of Prandtl 
number in the new model, we have one more free parameter to be adjusted. With the variation of 
this parameter, a continuum dynamic performance from ES-model to S-model, and beyond, can 
be identified. 



In the past years, a unified gas kinetic scheme (UGKS) [12l4l4i| has been well developed. The 
BGK model and the S-model have been employed in the UGKS. In this paper we will use the UGKS 
framework to construct numerical scheme for the generalized kinetic model. The numerical scheme 
will be used to exam physical performance of different kinetic models with the variation of the 
parameter, where both the ES-model and S-model become limiting cases. A continuous dynamic 



transition between these two models can be obtained. Through investigations, the performances 
of different kinetic models in the transition regime are presented in details. 

This paper is organized as following. Section 2 presents the UGKS for the ES-model and other 
kinetic models. Section 3 proposes a generalized kinetic model. Section 4 gives the simulation 
results of the new model in the shock structure and microti ow computations. The parameter 
dependent dynamic effect will be discussed in different test cases. Section 5 presents the analysis 
and insight of the new model. The last section is the conclusion. 

II. UNIFIED GAS KINETIC SCHEME FOR KINETIC MODELS 

The unified gas kinetic scheme is a direct modeling method to simulate flows in the whole 
Knudsen number regimes. It is a finite volume conservation law for the evolution of gas distribution 
function. Besides the evolution of conservative flow variables, such as density, momentum and 
energy, the time evolution of gas distribution function at discrete particle velocity is solved as 
well in order to capture the non-equilibrium molecular transport. Therefore, how to evaluate the 
fluxes of a gas distribution function across a cell interface is a central ingredient in UGKS. The 
kinetic model is always employed in UGKS to provide the evolution dynamics of the distribution 
function, but the UGKS is not targeting to solely solve the kinetic model itself, because the physical 
modeling scale of the kinetic model can be different from the numerical cell size scale. The UGKS 
is a direct physical modeling of flow motion in the scale of the discretized space and the integral 
solution used from the kinetic model covers the flow evolution from kinetic to the hydrodynamics 
scales. The specific fiux used at the cell interface depends on the ratio of time step to the local 
particle collision time. 

In this section, a brief review of the UGKS is presented. Since the Shakhov model has been im- 
plemented in UGKS, this section will introduce the UGKS with a general kinetic model. Generally, 
a kinetic model takes the following formulation. 

Of ^ df g+-f ,.- 

The / represents the velocity distribution function, and the g^ is the post collision term. The x 
and u represent the physical space variables and the velocity space variables respectively. Here r 
is relaxation time. 

The macroscopic quantities, such as, the mass p, momentum pJJ (pUi), energy pE, stress tensor 



P (pij) and heat flux q (qi), can be derived from the distribution function /, 



W 



pU = ipfdu, 
[pEj 

Pij = {ui- Ui){uj - Uj)fdu, 

qi = J^{u,-Ui){u-Uffdu, 



(2) 



where ip is defined as following, 



V' = (l,u,-u^) 



2\T 



(3) 



and du is the volume element in the velocity space. Since mass, momentum, and energy are 
conserved during particle collisions, / and g^ satisfy the conservation constraint. 



J{g+ - f)^du = 0, 



(4) 



at any location and any time. 

Taking the collision time as a local constant, there is an analytic solution from kinetic model, 

/(x,t,u,0=e-*/Vo(x-ut) 

+- r5+(x',t',u,0e-(*-*')/-dt', (5) 

T Jo 

where x' = x — u(t — t'). 

Applying this solution at cell interface, the mass flux, momentum flux and energy flux can be 

obtained as the following, 

T = T = 

macro ■J momentum 

\ •'energy I 

The ri„ denotes the entire velocity space. 

The flux of velocity distribution function at particle velocity u^ takes the following form: 






(6) 



Jufe = n • / u/du. 



(7) 



where 1^^^, denotes the velocity space around u^. The first term on the right hand side of Eq.Q 
can be directly evaluated from the initial distribution function. For simplicity, the cell interface is 



assumed to locate at x = 0. The normal direction of the cell interface is denoted by n. Suppose 
the initial distribution function takes the following form at a cell interface: 

fuf^ „ .X| _ fn (^^ _ 1 /o,fc(0) + "aF ■ ^' ^ - 0' fs.s 

where nonlinear limiter is used to reconstruct f^, f^, and the corresponding derivatives. 

The second term of Eq.([5]) corresponds to the hydrodynamic scale physics which should be 
constructed from macroscopic quantities. We can use a continuous distribution function to evaluate 
the integral term. For an equilibrium state g~^ around a cell interface, it can be formally expressed 
as following, 

5+(x,u,f)=g+ + r7+-x + 5+t. (9) 

In fact, the spatial and temporal derivatives of g~^ are the key components for the construction 
of UGKS. The derivatives of g~^ may be complicated. Fortunately, only in the continuum regime 
this term becomes important. Here, we use the derivatives of a local Maxwellian distribution to 
approximate these quantities. And we have, 

g+ (x, u, t) = g^ + 5(o,x • x + go^tt, (10) 

where g^ denotes the post collision state at the beginning of each time step at the cell interface 
and go is a local Maxwellian distribution function located at x = 0. It can be written as. 



vr 



K+2 



x\ — 



90 = P{- e-^^^-^r^ (11) 



where A = p/2p. 

For a specific kinetic model, g^ and go are uniquely determined by the initial distribution 
function f^. For example, the conservative variables are evaluated by applying the compatibility 
condition. The conservation constraint at (x = 0, t = 0+) gives 

Wo = J foH^ = Y,fo,ki^ 

= Y.(fo,kH[^ . Ufc] + /o^,(l - H[n . Uk]))^, (12) 



where H[x] is the Heaviside function defined by 

H[x] = 



0, X < 0, 

1, x>0. 



(13) 



Similarly, the high order moments, say, the stress tensor and the heat flux can be derived by Eq.([2]) 
at a cell interface. For the details of the numerical reconstruction, please refer to the articles about 



gas kinetic scheme [l^, |l3, UMj ■ 



Applying the conservation law, the evolution of flow quantities can be obtained. Owing to the 
absence of source term, the evolution of the macroscopic conservative quantities becomes, 

.*" + ! 



w 



n+1 



^" + T^/" yZ^SmJ'macrodt, (14) 



where V^- is the volume of f^x^ in the physical space, AS'm is the area of interface and m is index 
of surfaces of f^x^ • 

The collision term must be considered for the update of the distribution function. Here, we use 
two steps to update the distribution function. 

1 /■t"+^ n+(") _ fn 



m 



V-v /in 



+^( '"^ Jr + ^' . "^ )> (15) 



At flrst, we derive /* as a medium state. And then solving the second equation, we get f^^^ 
at the next time level. The above procedure is identical for an arbitrary g'^. For ES-model, g^ is 
written as 

g+ = g[f] = P exp(-i(u - U) • T-^ • (u - U)). (16) 

ydet(z7rij ^ 

Here, T is a tensor related to the stress tensor P, 

T = (1 - Ces)RTl + Ces^/p, (17) 

where R is gas constant and T is gas temperature. Andries provided a simple proof that ES-model 
preserves a correct Prandtl number [3]. The same proof can be done for Shakhov model. In the 
Shakhov model, the g^ takes the form, 

5+ = M[f]{l + (1 - Cshak)c ■ q(|^ - 5)/(5piiT)), (18) 

where A^[/] denotes the Maxwellian distribution function, T is temperature, q is heat flux, c = 
u — U is peculiar velocity and Cghak is a parameter which is related to the Prandtl number in this 
model. 



III. A GENERALIZED KINETIC MODEL 

Here we discuss the different ways to fix tlie Prandtl number in tlie kinetic models. Following 
Andries' proof [^l, we expand the distribution function in continuum regime, 

f = g+-r{M[f]t + u.M[f]^) + o{T). (19) 

Let's consider the ES-model first. As odd moments of peculiar velocity of Gaussian function is 
zero, the ES-model has no contribution to the heat flux of the distribution function. Therefore, the 
Prandtl number is effected only by the variation of stress tensor. The second term, —T{M[f]t + 
u • A^[/]x), corresponds to the contribution of BGK model to the stress tensor. The second order 
moments of peculiar velocity of Eq. (jl9p can be written as the following. 



P = {l-Ces)pRTI + CesP + 0{T)tgk- (20) 

Here, the definition of the stress tensor has been considered. And solving the P, we get, 

1 



1-a 



-0{T)bgk (21) 



Here, p is the shear stress defined as p = P — pRTI, and the 0{T)i,gk corresponds to the shear 
stress from derivative of local Maxwellian distribution function which is exactly the shear stress of 
the BGK model. The q from Eq. (119p will be identical to that in the BGK model. So the Prandtl 
number of ES-model is, 

Pr=:-4^Pr,,, = -^. (22) 

For the S-model, according to Eq. p^ , the heat flux of distribution function 

q = (1 - Cshak)(i + 0{T)bgk, (23) 

where the 0{T)hgk corresponds to the heat flux from the BGK model. And the q is 

1 



^ ^ -^bgk- (24) 

The Shakhov model does not affect the second order moments. So the stress tensor of Shakhov 
model keeps unchanged in comparison with the BGK model. So the Prandtl number for Shakhov 
model is, 

^'^shak = ^shak'^^'^bgk = ^shak- (25) 



These proofs imply that it is sufficient to achieve a correct continuum hmit as long as the spatial 
and temporal derivatives are expressed as the expansion of local Maxwellian. So, it is appropriate 
for the hydrodynamic flux to be estimated by the derivative of the Maxwellian function in the 
integral solution as mentioned in the last section. Furthermore, the proofs also show that the 
kinetic models fix the Prandtl number via the adjustment of either stress or heat flux of the 
relaxation term. It is quite straightforward to combine these two approaches together. 

The ES-model and S-model change either the stress tensor or the heat flux of the post collision 
terms to achieve a correct Prandtl number. How about to change these two quantities simultane- 
ously. It's obvious that this kind of modification could also generate a correct Prandtl number, 
and provides a free parameter as a by-product. 

Specifically, the post collision term of the generalized kinetic model is 

g+=g[f]+S[f], (26) 

where the g[f] is defined by Eq. (ll6p . The S[f] is Eq. (ll8p for Shakhov model without the first 
equilibrium state, namely, 

5[/] = Mimi - Cshak)c ■ q(-^ - 5)/{5pRT)]. (27) 

The two coefficients, Ces and Cshak, are two independent parameters at this moment. In order 
to obtain the right transport coefficients, we still follow the proof of Andries [3|. Eq. (|19p changes 
to the following one, 

/ = g[f] + S[f] - riMt + n-M^) + o{t). (28) 



For stress tensor. 



then 



P = {l-C,s)pRTI + C,sP + Pbgk, (29) 



Pbgk- (30) 



(1 - Ces) 



And for heat flux, 



then 



q = (1 - Cshak)<i + (Ibgk, (31) 



■^ Qbgfc- (32) 

^shak 
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As a result, the Prandtl number for the generahzed kinetic model is 



And the viscosity is 



Pr = p^Prftgfc = p^. (33) 



If Prandtl number is fixed, there is a free parameter in the generalized model. Here, the Ces can 
be taken as the free parameter. When Ces = and Cghak = Pr j the generalized model is identical 
with the Shakhov model. When Ces = 1 — -p- and Cshak = 1, it gives the ES-model. When 
Ces = and Cshak = 1, it presents the BGK model. And for the other values, the generalized 
kinetic model shows how the ES-model changes to Shakhov model continuously. And the new 
free parameter might provide an opportunity to preserve additional physical properties in the full 
Boltzmann collision term. 

The generalized kinetic model is employed in the UGKS introduced in section |TT] for its numerical 
solution. 

IV. NUMERICAL RESULTS 
A. Shock structure 

The shock structure is a typical example of non-equilibrium flow structure, and is a distinguish- 
able test case. Kinetic models show very different performances in shock structure simulation. To 
examine capabilities of ES-model and S-model, Mach 8 argon shock structure is simulated, and the 



solutions are compared with the DSMC results 



m- 



The DSMC code is provided by G.A. Bird. 



The viscosity-temperature coefficient a; is 0.81, namely, // ~ j^o.si^ -pj^g Prandtl number is 2/3, and 
Ces varies from —0.5 to 0.5 in the generalized kinetic model. The two special cases, the ES-model 
and the S-model, are included in this set of simulations. The reference viscosity is determined as 
following. 



30 pXV27rRT 

^'■^^"(7-2u;)(5-2u;) 4 ' ^^^' 

In our simulation, the spatial coordinate is normalized by the upstream mean free path, namely, 
the upstream mean free path of argon is 1. The computational spatial domain is [—50,30] and is 
uniform meshed by 300 grids. 
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FIG. 1. The shock structure for ES-model and S-model at Ma = 8 and w = 0.81. 
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FIG. 2. The shock structure for the generahzed kinetic model with different Ces at Ma = 8, and cj — 0.81. 



Figure [T] gives the density and temperature profiles from the ES-model and S-model. The S- 
model and the DSMC present almost identical density profiles. But the temperature rises a little bit 
early for the S-model. In comparison with S-model, the ES-model predicts a narrow density profile 
and a wide temperature profile. Obviously, the S-model performs much better than ES-model in 
this case. 

Figure [5] shows the tendency how the shock structure changes while the Ces varies from —0.5 to 
0.5. Note that Pr = 2/3 is fixed in all these results. The generalized kinetic model presents a set 
of shock structures with the same Prandtl number. When the Ces = —0.5, the generalized kinetic 
model presents the ES-model. As the value of Ces becomes larger, the temperature profile becomes 
steeper. Meanwhile, the density profile grows wider. When Ces is larger than 0, the temperature 
profile still becomes steepening. But when Ces exceeds 0.15, the density profile turns out to be 
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twisted near the upstream. Although the annoyed twisting density profile makes this range of the 
free parameter unacceptable, the strong dependence of the Cgs is confirmed. This coefficient effects 
the behavior of kinetic model. 

Taking moments of the generalized kinetic model, consider the following three equations for 
different moments, 

df 1 



a ;(/-9^). (36) 



dPjj ^ dpij_ ^ -(1 - Ces) 

dt ~ dt ~ T 

dqi —Cshak 



Pij, (37) 



at r *■ <="" 

which determine three different relaxation processes, namely, the relaxation of distribution function 
itself, the relaxation of second order moments and the relaxation of third order moments. The 
ratios between different relaxation rates are determined by the two coefficients, Ces and Cshak- 
Let the Prandtl number fixed, i.e., Pr = Cshak/{^ ~ Ces) keeps constant when changing Ces- The 
(1 — Ces) gives the ratio between the relaxation of distribution function and the relaxation of the 
second moments of distribution function. This is the physical meaning of the Ces- For example, if 
(1 — Ces) is bigger than 1, the second order moments decease more rapidly than the distribution 
function itself. As shown in figure [21 different Ces presents different relaxation ratio and provides 
different shock structures. 

B. Force driven Poiseuille flow 

In the force driven Poiseuille flow, the external force drives the flow motion between two fixed 
plates. The flow fleld will achieve a steady state when the external force is balanced by the shear 
stress from the flxed boundaries. We also consider monatomic gas in this simulation. To follow 
the study in [17|, the Knudsen number is deflned as, 



j^^^ /7/WM, (39) 

where L is the width of the channel, and subscript denotes the initial value of variable. The 
gas is conflned between two vertical plates which locate at a; = —0.5 and x = 0.5 respectively. 
The temperature of the plates is T^„ = 1. The initial flow states are shown as following, Tq = 1, 
Po = 1; Po = 1- The gravity is represented by G, and is in the vertical direction. Here the hard 
sphere molecule is adopted, namely, the viscosity-temperature coefficient is uj = 0.5. The gas- wall 
interaction uses fully diffusive kinetic boundary condition. Due to the large value of G = 1, this 
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FIG. 3. The velocity and temperature profiles from the generalized kinetic model under different Knudsen 
numbers. The Knudsen numbers are 1, 0.1 and 0.05 respectively from top to the bottom. And the Gravity 
isG = l 
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test can become a very tough one and the distribution function is fully distorted by the external 
forcing, especially at high Knudsen number. 

As shown in figure[3l when the Knudsen number is small, say, Kn = 0.05, the difference between 
results from different kinetic models and the DSMC is small. But, as the Knudsen number becomes 
large, the temperature profiles separate from each other. Similar to the shock structure, for all 
Knudsen number, the results from S-model are closer to the DSMC resluts than the ES-model. 
The profiles with different Ces cover the results of ES-model and S-model. And when the Ces is 
larger than 0, the temperature profile moves from the S-model result to the DSMC result. It is 
clearly shown that the generalized kinetic model can predict more accurate results in comparison 
with the S-model and ES-model if Ces is specified properly. 

C. Unsteady boundary heating 

In this section we solve the unsteady fiow problem using the generalized kinetic model. The 



numerical configuration is identical to the unsteady boundary heating problem in reference |l8l |. 
The gas is heated by two wall with time-dependent temperatures T^ = 1 + 0.002 sin(0t). Hard 
sphere molecule is adopted in the simulation. And the Prandtl number is 2/3. The Knudsen 
number is defined as following, 

5^/2^ PoL 



Figure [4] presents the velocity and temperature profiles at 9t = 3tt/2. The U velocity is normalized 
by 2 X 10~^, and AT is defined as AT = (T — To)/0.002. Obviously, the results from S-model is 
closer to the LVDSMC results. When C^s is larger than 0, the generalized kinetic model gives a 
better result. This coefficient is very close to the one in the force driven Poiseuille flow. 

D. Response of a gas to a spatially varying boundary temperature 



k 



varying boundary temperature 



19l |. Gas is confined between two 



The last simulation is about the response of a gas to a spatia 
in 2-D domain. The numerical setup is the same as the case in 
horizontal boundaries. The lower boundary at y = is fully diffusive with a temperature given by 
T^ = Tq{1 — 0.5 cos(27rx)). An identical boundary is located at y = 1. The Knudsen number based 
on the separation between the two boundaries is Kn = 1. Working gas is argon with reference 
viscosity defined by Eq. (135p . Owing to the symmetries in the x and y directions, the simulation 



domain is chosen as [0, 1/2] x [0, 1/2]. Figure [5] shows the temperature contour from ES-model and 
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FIG. 4. The velocity and temperature profile for unsteady boundary heating problem at Qt — "ii: j2. The 
is defined as = Tr\/2/16. 



S-model. The background data is extracted from reference [191 ]. Unlike the previous two test cases, 
ES-model predicts more accurate results in comparison with DSMC results. 
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FIG. 5. The temperature contour of a spatially boundary temperature variation problem. The dash line is 
DSMC data extracted from the reference [l9|. 

V. ANALYSIS 



The individual ES-model and S-niodel were constructed to set the Prandtl number as a free 
parameter. Physically, the Prandtl number of a monatomic gas has a fixed value, especially in the 
continuum flow regime. Therefore, the Prandtl number should not be taken as a free parameter 
for monatomic ideal gas. With a fixed Prandtl number, theoretically there is not any freedom in 
the ES-model and S-model. In this study we proposed a generalized kinetic model. Besides a fixed 
Prandtl number for monatomic gas, the new model provides one more free parameter. This free 
parameter can present a continuum spectrum of kinetic models with correct Prandtl number. This 
parameter provides ways to mimic more complicated physical relaxation process. With certain 
choices of this free parameter, say Ces, the S-model and ES-model become a subset of the new 
model. In the force driven Poiseuille flow and unsteady boundary heating problem, the new model 
provides a way to get accurate results when C^s is set to be larger than 0. 

As mentioned in last section, C^s and Cshak are related to the relaxation of moments of the 
distribution function. To shed light on this topic, we exam the Boltzmann collision term for VHS 
molecule. The Boltzmann equation is written as following, 

dinf) , dinf) 



+ u 



J{nf), 



(41) 



dt dx 

where / is normalized distribution function, n represents the particle number density, and J{nf) 
denotes the Boltzmann collision term. The collision integral is defined as 



A[Q] 



OO /■ + 00 l-i-K 



oo JO 



n'^Qiffi - ffi)cradndudui, 



(42) 
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where is the soHd angle for scattering molecule, c,. is the relative velocity between two colliding 
molecule, and a is the collision cross section. Consider a spatially homogenous monatomic gas 
problem. The moments equation of the Boltzmann equation gives the relaxation process of the 
moments. For the quantity Q, the relaxation process can be written as 

md<nf,Q> 

— = mA[Q], (43) 

where < nf, Q >= J nfQdu, and m is the mass of molecule. For example, if Q = u^ and for 
Maxwell molecule, i.e., /i ~ T, the corresponding relaxation equation is 



dPn dpii 21 

—- — = —- — = ml\\u . 



at m '-'-■■ (^) 



The collision integral can be obtained explicitly for Maxwell molecule [161], such as 

-^- = -Pii- (45) 

at II 

For other molecules, there is no explicit solution. However, some qualitative results can be deduced 
from a given distribution function. Here, we consider two kinds of distribution functions for VHS 
molecule. The diameter of VHS molecule is given by 

d = dref{Cr,ref/CrT, (46) 

where v = uj — 1/2. The first distribution function is the one employed in Grad's thirteen moments 



method 



2m\i and it reads 



/ = A^[/](l + (u-U).^.(u-U) 

+^-(u-U)(("-^)'-l)V (47) 

pRT ^ '^ 5RT ' j ^ ' 

For the case when (P — pI)/{2pRT) and q/i^pRT) are much less than 1, by substituting Eq. ([47ll 
into collision integral (Eq. (P2]) ). the above distribution function gives 

A[uu] (48) 

= -(n/m)^:^^-i^4-(i?r)^-r(4 - v)p, 
2 Iot/tt 

A[iuu2] (49) 

= -(n/m)^!^ J^4-(i?r)V-^r(4 - .)^q, 
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where F denotes the Gamma fmiction. The viscosity of the VHS molecule and the mean colli- 
sion rate (l/r) per molecule in an equilibrium gas of VHS molecules are given by [16l |. Here we 
reformulate them as following, 

Anc%^jaref'^~''{RTf'^-''T(2-v)/^ (51) 



r 

Using the above results, the relaxation process of moments of the Boltzamnn equation can be 
written as, 

^""^ --{nf-irJinf)+nf)), (52) 



dt T 

Pij, (53) 



dpij p 



dt /i 

ITT = -o-Qi- (54) 

dt 2, jjL 

Actually, for VHS molecule in a local equilibrium state, the Ces can be derived as |l6l |. 

1 30 p 



T (7-2a;)(5-2w)/x' 



(55) 



C. = 1-(1^M(^^M. (56) 

Here Ces is confined in a domain of [0.2, 0.5] for VHS molecules, and can be taken as a constant. 
However, there is no universal conclusion. For the shock structure calculation, the new model with 
such a range of Ces seems to give inappropriate solutions. 

Hereafter we consider another distribution function. Assume the distribution function is com- 
posed of two delta function, say, 

/ = a6{u — (1 — a)uQ) + (1 — a)5{u + auo), (57) 

where a G [0,1], u denotes molecule velocity in x direction. Then the pressure, stress tensor and 
heat flux can be expressed by a and uq, 

p = -mna{\ — a)u'^, (58) 

Pii = 3p, (59) 

qi = -mna{\ — a)(l — 2a)u^. (60) 



The collision rate is 



- =2n^d2^^c2;;^^a(l - a)nj~2'^, < t; < 1/2, (61) 

- = n-Kdl^jCr^ref, i; = 1/2. (62) 
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Note that the 1/r is not continuous when v = 1/2, because two molecules with identical velocity 
collide with each other with infinite collision cross section between them. It is inappropriate to 
count this kind of collision. So we will discuss the case of < f < 1/2. Substituting Eq. (f57ll into 
Eq. (|42p . the collision terms give, 



mA[u ] 

mA\—uu 1 

^2 ^ 



1 (Qi/P? + {Pu/pf 
-Pn- 



T 



{Pll/Pf 



\2^ {qi/pf + {Pll/pf 



r 3 



(63) 
(64) 



{Pii/P? 

Here, the relaxation process is totally different from the near equilibrium state as shown before. 
Ces in this case can be formally written as 



a 



1 



{Pii/p? 

Obviously, it is not a constant. Furthermore, it is less than and can even go to minus infinity. 
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FIG. 6. The shock structure from the generalized kinetic model with a variable Ces at Ma — 8, and ui = 0.81. 



With this understanding, we construct a variable Ces in the shock structure calculation in 
order to get a good agreement with DSMC. As show in figure El a perfect shock structure can be 
obtained and the corresponding Ces is plotted for this calculation. The temperature profile is much 
improved, while the density profile changes only a little bit. The early raising of temperature in 
the upstream is suppressed efficiently. 

For boundary temperature variation problem, the value of Ces is preferred to recover the ES- 
model, namely, Ces = —0.5. The figure [7] shows the distribution function at (0,0.5). Based on 
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the above analysis, two peak structure corresponds to a negative value of Ces- Therefore, we 
qualitatively conclude that the ES-model is more appropriate for this problem. Based on these 
numerical results and analysis, we believe that this new free parameter has significant physical 
insight which deserves its further study. 




FIG. 7. The distribution function at location (0,0.5) for spatially varying boundary temperature heating 
problem. 

As mentioned above, the relaxation rate of different moments depends on the distribution 
function and molecular types. And Ces cannot be taken as a constant for transition flow. For 
the two-peak distribution functions, this coefficient could be even far less than 0. But Ces in the 
ES-model is always constrained in the interval [—0.5, 1) in order to keep a positive eigenvalue of 
T. In fact, an alternative of Gaussian distribution can be adopted in the kinetic model, 

g[f] « M[m + (u - U) • r • (u - U)), (65) 

where 



I 

2{Krf 



^i = ?77FFW^^J' ^^^'' (^^) 



and 



Surprisingly, although the above expansion cannot guarantee the positivity of the distribution 
function, the numerical results from the above expansion are very close to that where a full Gaussian 
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distribution function is used. It indicates that for any formulation we adopt in the model equation, 
the results of macroscopic variables will be the same, as long as the moments of the collision term 
are identical from different kinetic models. Furthermore, replaced by the expansion, the lower 
bound of Ces for Gaussian distribution can be removed. We can use a value of Ces less than —0.5. 

VI. CONCLUSION 

In this paper, we have developed a generalized kinetic model through the combination of ES- 
model and S-model. With a fixed Prandtl number, this new model provides an additional free 
parameter, which can be used to recover the physical solution more accurately. By changing the 
free parameter the different relaxation time between different moments of a distribution function 
can be simulated. The unified gas kinetic scheme is used for the construction of numerical solution 
of the generalized kinetic model. With the variation of this free parameter, the new model covers 
the BGK model, ES-model and Shakhov model. At the same time, it provides a continuum 
spectrum of kinetic models and different dynamics with a variation of this parameter. In most 
cases, the S-model presents more accurate numerical results. The numerical study indicates that the 
essential property for a kinetic model to capture physically valid solutions is the ratios between the 
relaxation rates of different moments of a distribution function. We believe that the introduction 
of this generalized kinetic model is important in the study of non-equilibrium flow, and this free 
parameter has significant physics basis, which deserves its further study. 
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